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Abstract 

We present unquenched lattice QCD results for the matrix elements of four-fermion operators 
relevant to the description of the neutral K and D mixing in the Standard Model and its 
extensions. We have employed simulations with Ny = dynamical sea quarks at three 

values of the lattice spacings in the interval 0.06 - 0.09 fm and pseudoscalar meson masses 
in the range 210 - 450 MeV. Our results are extrapolated to the continuum limit and to the 
physical pion mass. Renormalization constants have been determined non-perturbatively in 
the RI-MOM scheme. In particular, for the Kaon bag-parameter, which is relevant for the 
mixing in the Standard Model, we obtain = 0.717(24). 




1 Introduction 


Lacking experimental data above production threshold, flavor physics offers the unique 
possibility for an indirect discovery of New Physics (NP) effects through virtual exchanges 
of yet-to-be-discovered heavy particles in loop suppressed processes. This approach, which 
is particularly promising for processes that are highly suppressed within the Standard Model 
(SM), proved to be very successful in the past allowing for the indirect determination of the 
charm and top quark mass |[T]-j^. 

Moreover flavor physics data play a major role in providing stringent tests of the CKM 
paradigm and allowing the determination of the magnitude of the mixing matrix elements. 
In particular AS = 2 and /S.B = 2 flavor changing neutral current processes are crucial to 
the unitarity triangle analysis. They are also quite valuable in constraining NP models, see 
e.g. i|5{n} with data on AS = 2 oscillations providing the most stringent constraints m 

m- 


Of special interest are the AC = 2 transitions occurring in oscillations 113 - 


and |[T^[n|, as this is the only SM process in which mixing involves up-type quarks. 
CP violation through these mixings is expected to be strongly suppressed within the SM, 
because they are dominated by light (d, s) quark exchange entailing also important long 
range interactions. Thus any experimental signal of CP violation in the neutral D meson 


sector would be a strong indication for the existence of NP 118 -221. Even in the absence 
of CP-violation, our determination of AC = 2 operator matrix elements allows to put 
constraints on Beyond the Standard Model (BSM) models. 

In this paper we present a determination of the bag-parameters relevant for the de¬ 
scription of the AS* = 2 and AC = 2 transitions. We compute meson-anti-meson matrix 
elements of the whole basis of dimension-six four-fermion operators contributing the most 
general form of the effective AF = 2 Hamiltonian |[23]-(26|. Beyond the “left-left” oper¬ 


ator, relevant for the SM, flavor-changing extra terms appear. The full effective AF = 2 
Hamiltonian reads 




i=l 


i=l 


( 1 ) 


where C* and Ci are the Wilson coefficients and the bare operators, Oi and Oi, correspond¬ 
ing to renormalized operators appearing in Eq. ([T]l, are 

Oi = - 75)^“] - 75)^^] , Oi = [d" 7 ^(l -h 75)^“] [^^ 7/^(1 + 75)^^] , 

O2 = [^“(1 - 75)^“] [^^(1 - 75)^''] , 4 = [^“(1 + 75)^“] [^^(1 + 75)^''] , 

O3 = [^“(1 - 75)^^] [h^{l - 75)^“] , O3 = [^“(1 + 75)^^] [h^{l + 75)^“] , 

04= [/i“( 1 - 75 K 1 [d^(l+ 75 K^] , 

05= [^“( 1 - 75 K^] [V(1 + 75 K 1 , 

(2) 

where a, /3 are color indices. The operators Oi_3 are obtained from Oi_3 with the replace¬ 
ment of (1 — 75) —)• (1 -|- 75). Since Oi_3 and Oi_3 have identical parity conserving 
parts, parity invariance of QCD allows to restrict our attention to only the set of oper¬ 
ators Oi, i = 1 ,..., 5 . In the present work we focus on the cases (/i, i) = (s, d) and 
(hj) = (c,u). 

The Wilson coefficients describe short distance effects. Accordingly, they will also 
depend on the heavy degrees of freedom possibly circulating in loops. The low energy 
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dynamics is incorporated in the matrix element of the operators Oi. The renormalization 
scale ^ gets compensated between the Wilson coefficients and the matrix elements of the 
renormalized operators. 

Lattice QCD provides a first principles determination of the bag-parameters Bi. These 
are dimensionless quantities defined as fhe ratio of fhe non-perfurbafively compufed four- 
fermion mafrix elemenf over fhe value this matrix element takes in the vacuum saturation 
approximation. The reason for working with ratios is that they offer the advantage of a 
substantial cancellation of systematic and statistical uncertainties between the numerator 
and the denominator. For their definition (see for example Ref. p7| ) 

(pVi(/i)|pO) =6i?i(F)m^o/|-o 

for f = 2,...,5, 

+ mh{n)) 

where = {8/3, —5/3,1/3, 2, 2/3} and stands for either a or a pseudoscalar 
meson. The corresponding mass and decay constant are denoted bjj^mpo and fpo, respec¬ 
tively. The quantities and mh{ii) are the light and heavy quark masses of the neutral 

K and D pseudoscalar mesons are made of, renormalized at the scale /r. 

We see from Eq. ^ that, while the matrix element {P^\Oi{^)\P^) vanishes as the 
pseudoscalar mass goes to zero, the four other matrix elements do not. We recall that Bi 
parametrises the SM operator while the Bi, i > 2, parametrise the BSM ones. 

The computations presented in this paper have been performed making use of the Nf = 

2 -I- 1 -|- 1 dynamical quark gauge configurations generated by the European Twisted Mass 
Collaboration (ETMC) | [28][29| at three values of the lattice spacing, a ~ 0.06 — 0.09 fm, 
with the lightest pseudoscalar mass values in the range Mpg ~ 210 — 450 MeV. Spatial 
lattice sizes areL~2.1 — 3.0fm with MpgL ~ 3.1 — 4.5. Operator renormalization has 
been performed non-perturbatively in the RI-MOM scheme . 

Results - Eor the reader’s convenience we immediately summarize our main results 
for and meson bag-parameters. 

We collect in Table [T] the values of the five bag-paramefers fhaf are required fo describe fhe 
neufral Kaon mixing in fhe SM and beyond. We give fhe numbers in fhe MS renormalizafion 
scheme of Ref. | |^ and in the RT scheme at the scale of ^ = 3 GeV. Eor results given in the 
MS scheme the second quoted error provides our estimates for the systematic uncertainty 
coming from the perturbative matching between RT and MS schemes, which range from 
0.5% to about 4%. The uncertainties on the central values stemming exclusively from our 
lattice computations are given by the first error and range from about 3.3% to 7.5%. 

Neglecting the tiny over-unquenching error due to the presence of the charm in the 
sea (see discussion below) we adopt the continuum limit results in the MS scheme at /r = 

3 GeV as our best estimate of the desired Kaon mixing bag parameters in QCD with u, d 
and s active flavors for the same scheme and scale. Then for Bk = Bi we find in the RGI 
scheme the value R 

^RGi{Nf=3) ^ Q 7^7(24), (4) 

with a total uncertainty of about 3.4%. 

'in our formulae we use the notation for a neutral pseudoscalar meson although we work in the isospin 
symmetric limit; so in practice we make no distinction between the masses and decay constants of the neutral 
and charged pseudoscalar mesons. 

^If for converting our continuum limit Bi in the MS scheme at /r = 3 GeV to its RGI counterpart we had 
taken Nf = 4 we would have obtained = 0.728(24). 
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K°-K^ 

MS (3 GeV) 

0.506(17)(3) 

0.46(3)(1) 0.79(5)(1) 

0.78(4)(3) 

0.49(4)(1) 

RT (3 GeV) 

0.498(16) 

0.62(3) 1.10(7) 

0.98(5) 

0.66(5) 


Table 1: Continuum limit results for the bag-parameters Bj (i = 1,..., 5) relevant to the 


mixing renormalized in the MS scheme of Ref. 1311 and in the RT scheme at 
the scale of /x = 3 GeV. For results given in the MS scheme the second error indicates an 
estimate for the systematic uncertainty owing to the perturbative matching of RI' and MS 
schemes. 


In Table 1^ we summarize the results for the bag-parameters relevant for the case of the 
oscillations. For results given in the MS scheme the second error we quote rep¬ 
resents our estimate of the systematic uncertainty coming from the perturbative matching 
between RT and MS schemes. The uncertainties stemming only from our lattice computa¬ 
tions are given by the first quoted error and range from about 4% to 8%. 


d” - 

MS (3 GeV) 

0.757(27)(4) 

0.65(3)(2) 0.96(8)(2) 

0.91(5)(4) 

0.97(7)(1) 

RT (3 GeV) 

0.744(27) 

0.87(5) 1.34(11) 

1.14(6) 

1.39(9) 


Table 2: Same as in Tablej^for the mixing. 


The results of this paper are compared with the existing unquenched determinations 
inFigs.[T}|^ an d0 

In Fig. [T] we present a compilation of recent (unquenched) RGI values of Bk- Lattice 
computations are quite accurate with a total uncertainty of only a few percent. It is worth 
noting the rather weak dependence of B^^^ on the number of dynamical flavors. Our 
current B^^^ value compares well with = 2 -|- 1 and Nf = 2 lattice computations. 

Our result is also found to be in agreement with the estimate obtained using a model based 
on the dual representation of QCD as a theory of weakly interacting mesons for large N, 
which predicts negative sign corrections to the large N limit estimate given by the value 
Bk = 0.75 Il34)|35|. _ 

A comparison of recent determinations of the K —K^ bag-parameters Bi,i = 2,... ,5 
is presented in Fig. ^ . The ETM, RBC/UKQCD and SWME collaborations give for B 2 


and i ?3 results that are well compatible within the errors. A tension of up to 3 standard 
deviations is visible, instead, in the case of B 4 and B^ after the updated (preliminary) work 
of SWME 1^. 

Einally, in Eig. we show the comparison of the available results for the bag- 
parameters coming from ETMC computations with Nf = 2 and Nf = 2-|-l-|-l gauge 
configurations!^ 

Comparing the results for the bag-parameters collected in the Tables and with the 
ETMC results published in Refs 112| and |50|, the latter obtained with Nj = 2 dynami- 


^For recent reviews see Ref |32||33). _ 

'’For older quenched computations of the BSM Bi see Refs. j42||43). _ 

^Work in progress of an unquenched Nf = 2-1-1 computation for the D-mixing is reported in Ref. |47| 
For older works using quenched simulations see Refs j48[|49[. 
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Figure 1: A compilation of unquenched lattice results for the RGI value of the Bk pa¬ 
rameter. From top to bottom results are taken from Refs. |12|, |36|, |37|, |38|, |39| 
132|, |40|, |41|. Circle, squares and triangle correspond to Nf = 2, Nj = 2 + 1 and 
Nf = 2 + 1 + 1 dynamical quark computations, respectively. The full blue square indicates 
the FLAG average |32| over Nf = 2 1 data. For reader’s convenience some information 

on the basic features of each computation is also given. 
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Figure 2: A compilation of meson bag-parameters Bi, i = 2,..., 5. From top to bottom 
data have been taken from Refs p^|44[|45|. The work reported in Ref. 1441 is an updated 


computation of Ref. 1391. The label “CL” stands for continuum limit computation. Work in 
progress by the RBC/UKQCD has been reported in |461. 
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Figure 3: Results for the meson bag-parameters Bi, z = 1,..., 5 obtained by the ETMC 
with A^y = 2-|-1-|-1 (this paper) and = 2 |50| dynamical flavor lattice simulations. 


cal quark simulations, we notice that they are all compatible among themselves and have 
similar total uncertainties. Therefore, ceteris paribus, the main conclusions presented in 
these works concerning model-independent constraints on the NP scale from AS" = 2 and 
AC = 2 operators within the unitarity triangle analysis, remain unchanged. 

Before concluding this discussion, we find it useful to comment further about the de¬ 
pendence of Bk and of the other i?-parameters for mixing on the number of 

dynamical quarks. As known, the standard theoretical formula that provides the indirect 
CP violation parameter ck is obtained, through the low-energy effective weak Hamiltonian, 
after integrating out the heavy degrees of freedom including the charm quark. The reason is 
that it is only when the charm quark is integrated out, i.e. at scales p, ~ rUc, that the imag¬ 
inary part of the effective Hamiltonian for mixing becomes local (at the leading 

order in the 1/mc expansion). The advantage of this approach is that the long distance con¬ 
tributions to the amplitude, being related to the matrix elements of local operators, are more 
easily accessible to lattice computations. The price to pay, however, is that perturbation 
theory is uncertain at scales around the charm mass and, in addition, subleading corrections 
proportional to powers of p\/m^ (where pK = 0 {mK, ^qcd)) may not be negligible, 
particularly when aiming at a theoretical prediction for ck with percent precision. 

In the standard approach, both short-distance Wilson coefficient and long-distance ma¬ 
trix elements of the effective Hamiltonian have to be computed in the presence of three ac¬ 
tive quarks. In this respect, therefore, the lattice computation of H-parameters for 
mixing presented in this paper, being based on simulations performed with Nf = 2 -|- 1 -|- 1 
dynamical quarks, introduces a systematic error. Previous experience with Aj = 2-|-1-|-1 
lattice calculations suggests that the effect of the dynamical charm quark is presumably tiny, 
so that its impact in the determination of physical observables, which is undesired in this 
particular case, is likely too small to be detected at the level of the current precision. 

It should be also noted that a similar source of systematic error is introduced in the 
lattice calculations of Bk performed with only Nf = 2 -|- 1 dynamical quarks. In the 
latter case, indeed, the effect of the dynamical charm, which properly is not introduced 
in the determination of the matrix elements, is missing however in the lattice computation 
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of the hadronic observables which are needed to fix the action parameters. Therefore, the 
determination of the lattice scale as well as of the strange and light quark masses within any 
Nf = 2 + 1 lattice calculation is affected by the systematic error due to the quenching of the 
charm quark. This error then propagates into the calculation of Bk- These effects, namely 
the latter one and the error introduced in the A^y = 2 + 1 + 1 calculations of Bk, have the 
same physical origin and are presumably comparable in size. As already noted, both effects 
are likely to be currently negligible, as also indicated, a posteriori, by the good consistency 
observed in Fig. 1 among the lattice determinations of Bk based on different numbers of 
active quarks. 

In order to bypass this error, the theoretical determination of ex should be performed 
by keeping an active charm in the calculation. With the advances of the lattice technique, 
the computation of the matrix elements of non-local operators has becoming feasible and a 
first, exploratory lattice calculation of the real part of the effective Hamiltonian for 
mixing has been presented in Refs 151 521. The same technique can be also applied to the 
calculation of the imaginary part of the Hamiltonian, which is relevant for ex- While these 
lattice studies are not yet as accurate as the standard computations of local operator matrix 
elements, they are opening a new perspective and are likely to allow, in a near future, a 
significant improvement in the accuracy of the theoretical predictions of ex, for both the 
long distance and short distance contributions. 

Plan of the paper - In Section we review the simulation details and discuss our 
computational and analysis setup. Final results and a full account of the error budget is given 
in Section Finally, in the Appendices and we discuss the procedure we employed 
to determine in RI-MOM the full 5x5 renormalization constants (RCs) matrix for the 
four-fermion operators. 


2 Computational details 

In this work we have employed the mixed action twisted mass/Osterwalder-Seiler setup 
proposed in Ref. which provides automatic 0(a)-improvement and continuum-like 
renormalization pattern for the four-fermion operators, with only O(a^) unitarity violations. 


2.1 Lattice setup 


For the action of the light mass-degenerate sea quark doublet we have used the expression 
of Ref. 15^ which reads 


Sp = a 


X] I - “'•'s’’’ 




+ Msea ? 'Ipiix), 


(5) 

where it is intended that the untwisted mass has been tuned to its critical value, Mcr- As 
usual, the symbols and V* represent the nearest neighbour forward and backward co¬ 
variant derivatives, we define fhe quark doublef ^|;£ = '4’d)'^ and ^sea is the (light) sea 

twisted quark mass. 

With similar notations we take the action for the strange and charm quark doublet in 
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the sea | [55l to be 


S'/I = 2_/ ~ 


Mcr- - 




+ /i<7 + ^J-5T^j tphix), 

( 6 ) 


where /iff and are the bare twisted mass parameters from which the renormalized strange 
and charm masses can be derived. Pauli matrices in Eqs (|^ and Q act in flavor space. For 


more details on the twisted mass setup see Refs. |28 


Valence quarks are introduced via Osterwalder-Seiler (OS) fermions |591. The valence 
action is written as the sum of the individual quark flavor contributions in the form 


S°^ = a 


E 


E «/ 

f=i,l',h,h' 


1 


7/. (V^ + V*) - ij5rf 


M„- 




+ fif}qf{x). 


(7) 

where the label / is let to run over the different valence flavors f = £,£',h, h'. In the neutral 
K case light {£) and heavy {h) flavors denote down and strange quarks, respectively, while 
in the neutral D-case they stand for up and charm quarks. With the choice Vh = = r^,/ = 

—r^j one can prove |[5^ that at maximal twist automatic 0(a) improvement and absence 


of wrong chiral mixings 1601 is guaranteed. Flavor by flavor bare valence and sea quark 
masses are set equal to each other which is enough to keep unitarity violations to O(a^). 
The multiplicative mass renormalization constant i& Zp for all fermions. 

The lattice setup described above has been already successfully applied to determine 
the full set of four-fermion operator matrix elements relevant for the — K^, 

-f;0 


and oscillations in Refs. (1 2 


50 


61 


63 


2.2 Simulation Details 

We have used Nf = 2 + 1 + 1 gauge configuration ensembles, produced with the Iwasaki 


gluon action |64| and maximally twisted Wilson fermions, generated by the ETM Collabo¬ 
ration | |2^[29l |. 

In Table [3] we summarise the main simulation details relevant for the sea and valence 
sector. Simulation data have been taken at three values of the lattice spacing, namely a = 
0.0885(36), 0.0815(30) and 0.0619(18) fm, corresponding to /3 = 1.90, 1.95 and 2.10, 
respectively (see Ref. |[65|). 

As we said, light valence and sea quark masses are set equal, leading to pion masses in 
the range between 210 and 450 MeV. Strange and charm sea quark masses are chosen close 
to their physical value. To allow for a smooth interpolation to the physical values of the 
strange and charm quark mass, we have inverted the heavy valence Dirac matrix for three 
values, //“s”, of the strange quark mass and three values, /i“c”, of the charm mass, around 
the corresponding physical mass values. 

The lattice scale has been fixed using /,r- The u/d, strange and charm quark masses 
have been determined comparing with the experimental values of the pion, K and 
meson mass, respectively. Further details of our simulation setup can be found in Ref. |651. 

Valence light and strange quark propagators have been computed by employing spatial 
stochastic sources at a randomly chosen time-slice, adopting the “one-end” trick stochastic 
method of Ref. | [66p7| . In correlators where the charm quark is involved, Gaussian smeared 
interpolating quark fields [ [68| are used in order to suppress the contribution of excited states. 
This allows ground state identification at precocious Euclidean time separations. 
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For the values of the smearing parameters we take /cg = 4 and Nq = 30 Gaussian. 
In addition, we apply APE-smearing to the gauge links 1691 in the interpolating fields with 
parameters oape = 0.5 and Nape = 20. 


p 

P xT 

— ^/^sea 

aH“s" 

aH“c" 

1.90 (a-^ ~ 2.19 GeV) 

2P X 48 

0.0040 

0.0145 0.0185 0.0225 

0.21256 0.25 0.29404 

Her = 0.15 Hs = 0.19 


0.0060 





0.0080 





0.0100 




32=* X 64 

0.0030 

0.0145 0.0185 0.0225 

0.21256 0.25 0.29404 



0.0040 





0.0050 



1.95 (a-i ~ 2.50 GeV) 

24=* X 48 

0.0085 

0.0141 0.0180 0.0219 

0.18705 0.22 0.25875 

Ha = 0.135 Hs = 0.17 

32=* X 64 

0.0025 

0.0141 0.0180 0.0219 

0.18705 0.22 0.25875 



0.0035 





0.0055 





0.0075 



2.10 (a-i ~ 3.23 GeV) 

48=* X 96 

0.0015 

0.0118 0.0151 0.0184 

0.14454 0.17 0.19995 

Ha = 0.12 HS = 0.1385 


0.0020 





0.0030 




Table 3: Details of the simulation setup. Sea and valence fermion actions are displayed in 
Eqs. Q, Q and Q. 


2.3 Lattice operators and bag-parameters 


A detailed account of the lattice operators entering two- and three-point correlation func¬ 


tions was presented in Appendix [A| of Ref. |62|. Eor the reader’s convenience and to fix 


the notation we recall here some basic information. In our mixed action setup one needs to 


consider the following set of four-fermion operators 


0^4 = + [^7M759F][9f/7M759f/]) ± ^ C)] 

02[±\ = ± o t')} 

Ostti = ^ ^0} 


= 2 {([^ 9 “][ 49 £'] - ± 4 ^ /)} 


O5Y4 = ± ^ ^ 0 } > 


( 8 ) 


where a and j3 are color indices, the square parentheses denote spin covariant operator 
factors and the label “MA” stands for “Mixed Action”. 

We have set periodic boundary conditions for all fields, excepf for fhe quark fields 
which obey anfi-periodic boundary condifions in the time direction. Two “wall” operators 
with P°-meson quantum numbers (recall can be either or D^) are introduced at time 











slices yo and yo + Tsep/2. The first operator is constructed in terms of qi and quark fields 
and the second in terms of q^i and quark fields. Explicitly they are given by 


_ I ^ 

' yo ' - 


V, 


rh’ 

yo+Tsi 


L 

aN3 

L 


'^qe{y,yo)ibqh{v,vo) 

y 

X] ’ 2^0 + TsephbQh’ {y, yo + T, 


sepj 


(9) 


y 


In terms of them, the correlation functions we need to calculate are 

Ci{xo) = (I) +T.ep Oi|+f (x, xo) i = 

X 

Cxpixo) = (|)“ 2 (X"(I,W)P“), 

X 

C'xpixo) = (;^) ^C^yo+Tsep ^ (®! ^o)) 


( 10 ) 

( 11 ) 

( 12 ) 


where X can be either the axial current, Aq, or the pseudoscalar density, P. 

For three-point correlation functions with heavy (charm or heavier) quarks we can 
achieve reduced statistical uncertainties by decreasing the time separation between the two 


^). Therefore, while in the 


sources and using smearing techniques (see Refs. 
case we have set Tgep = T/2, in the case we have instead produced correlation 

functions setting Pep/a = 18 at = 1.9, T^ep/a = 20 at /3 = 1.95 and Psp/a = 26 at 
/3 = 2.10. 

Estimators for the bare bag-parameters are extracted from the asymptotic time be¬ 
haviour of the ratios of the three- to two-point correlators 


n^i\xo) = 


Ci(xo) 


C'ap(xo)C'^p(xo)’ 


^1=2 


>) ( \ — Ci=2,...,5ixo) 

- 2,...,5 0 Cpp{xo)Cpp{xo) 

(13) 

which for large time separations, yo 'C xq <C yo + Tsep, tend to the desired (bare) bag- 
parameters 


yQ<.xo<.yo+Tsep 

n[^\xo) -^ 




(p°Kio) (oiAri^°) 


(b) 


= 


(14) 


T^f\xo) 


ya^^o^yo+Tsep 
-)> 


/-FsO, 






(b) 




t = 2,...,5 (15) 


Figs 1^ and refer to the neutral K and D meson case, respectively. They illustrate 
the quality of the plateaux from which the estimates for the bare Bi {i = 1,..., 5) bag- 
parameters are extracted. The three panels correspond to three values of the lattice spacing 
at which simulations were performed. 

We note that our plateau choices are rather conservative and by reasonably varying the 
plateau interval (e.g. considering a plateau’s length as long as the double of the size of our 
principal choice) we find that the maximal systematic uncertainty is at the sub-percent level 
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(with maximal estimates being 0.2% for B 2 and 0.5% for for the neutral kaon and D 
cases, respectively). We therefore conclude that this systematic uncertainty is so smaller 
than the statistical one, indicated by the label “stat+fit+RCs” in Tables 4 and 5 (see below 
for details) that it can be safely neglected. 


(a) (I,) (b) 

S , 1.80 
P? l.GO 
1.20 
4 ? 1.00 

„ 0.30 

"pi 0.25 
_ 2.00 
1.90 
0.70 
0.65 
0.60 

^o/Tspp a^o/7sep 

Figured: Tl\^\xo) {i = 1,...,5) simulation points plotted against xo/Tsep for the 
case. Panel (a) (5 = 1.90, {a^i^a^s) = (0.0030,0.0185), volume = 32^ x 64. Panel 
(b) f3 = 1.95, {ane,ans) = (0.0025,0.0180), volume = 32^ x 64 . Panel (c) /3 = 2.10, 
{afi£, afis) = (0.0015,0.0.0151), volume = 48^ x 96. The dotted lines delimit the plateau 
region. Points for TZ^i \ ... at /3 = 2.10 have been slightly shifted upward by +0.05 
for accommodating data from all three /3’s in the same plotting scale. 
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Figure 5: Tlf\xo) {i = 1,..., 5) simulation points plotted against xo/Tgep for the — 
case. Panel (a) (3 = 1.90, (a^£,a^c) = (0.0030,0.25), volume = 32^ x 64. Panel 
(b) l3 = 1.95, (ane^afic) = (0.0025,0.22), volume = 32^ x 64. Panel (c) /3 = 2.10, 
{afi£,anc) = (0.0015,0.17), volume = 48^ x 96. The dotted lines delimit the plateau 
region. Points for and at /3 = 2.10 have been slightly shifted upward by +0.1 for 
accommodating data from the three /3’s in the same plotting scale. 
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2.4 Computation of the renormalized bag-parameters 


The renormalization pattern of the bag-parameters in our mixed action setup has been dis¬ 
cussed in detail in Refs. 112 ^62|. The renormalized bag-parameters are given by 
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Bi = 


^11 pjb) 

^iZaZv ^ 


B^ = 


p(fe) 

iiZpZs ^ 


hj = 2, ..,5, 


(16) 


where (^ 1 ,^ 2 , ■^ 3 , ^ 4 , '^ 5 ) = (8/3,—5/3, 1/3, 2, 2/3). The RCs of bilinear operators, 
namely Zy, Za, Zp and Zs, have been determined non-perturbatively in the RI'-MOM 


scheme in Ref. 1651. The four-fermion RCs, Z*, , have also been computed non-perturbatively 


in the same scheme. The calculation is presented in the Appendices [A] and [B| 

At each value of the light quark mass m = //sea our estimates of the bag-parameters 
are linearly interpolated to the physical strange (for neutral iT-mixing) or charm (for neutral 
O-mixing) quark mass. In both cases the interpolation turns out to be very smooth. A 
simultaneous chiral and lattice spacing extrapolation to the physical value of the pion mass 
and the continuum limit is finally performed. The u/d, strange and charm quark masses 
have been evaluated in the continuum limit in Ref. | |65t . 

Both for neutral K and D meson mixing studies and for all the Si’s we have employed 
a linear fit ansatz of the general form 


Bi — Bf -|- hijl^ -|- DiO^, 


(17) 


which in all cases nicely fits the data. In our notation the hat ('') symbol denotes renormal¬ 
ization in the MS scheme at the 3 GeV scale. We have also considered fit ansatze based on 


NLO ChPT 1701 for the K bag-parameters, of the kind 


S, = S 


1 , r - ‘^Bofli 2So/i£ 

1 -h bifie =F log 


16vr2/2 167r2/2 

and for the D bag-parameters the NLO HMChPT fit ansatze | |7T| 


DiU 


(18) 


Si = Sf 


B^ = Bf 


^ - (1 - 3^2) 2Bofie 2SoA£ 

1 -h ---log 


16vr2/2 167r2/2 


. , r ^ (1 T rg'^y) 2SoAr , 2SoAr 

1 + =F- 2 - mWi 


+ Dior 


+ Dia^ f = 2 ,4,5 . 


(19) 


In Eq. (18 1 the sign in front of the logarithm is minus for z = 1,2,3 and plus for z = 4, 5, 
whereas in the second of the Eqs. ( fTO] ) the sign is minus for z = 2 and plus for z = 4, 5. 
In the fit procedure we use the determinations of Sq and /o reported in Ref. . We also 
make use of the value Y = 1 derived in Ref. 0 and of the estimate g = 0.53(4) obtained 
from the lattice measurement of the gD*Dn coupling | |72| . In HQET the bag parameter 
S 3 is related to Si and S 2 . In particular, by setting Y = 1, S 3 and S 2 acquire identical 
logarithmic terms. 

In Eigs. and [7] we display the combined chiral and continuum fit of the Sj (z = 
1,..., 5) neutral K and D meson bag-parameters renormalized in the MS scheme of Ref. 1311 
at the scale of 3 GeV. We report in left and right panels, respectively, data corresponding 
to the two ways, namely Ml and M2, of determining the RCs proposed in Ref. | |73| , dif¬ 
fering in the manner 0 ( 0 ^) lattice artefacts are treated. More details will be given in the 
Appendices and Eor each Sj the results using either the polynomial or the chiral fit 
ansatze defined in Eqs. ( [TT] ), ( [T^ and ( [T^ are compafible among fhemselves wifhin less 
fhan one sfandard deviation. We also nofice fhaf fhe use of M1 and M2-lype RCs, fhough 
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Figure 6: Combined chiral and continuum extrapolation for the five Bi of the 
case. Bag-parameters are renormalized in the MS scheme of |311 at the scale of 3 GeV. 
Left and right panels correspond to Ml-type and M2-type four-fermion RCs, respectively, 
following the nomenclature of Ref. [ [73| . In each panel open circles and stars represent the 
value at the physical point corresponding to the linear and NLO ChPT fit, respectively. 


leading to rather different discretisation effects, provide compatible continuum limit deter¬ 
minations for the bag-parameters within 1-2 standard deviations. Moreover, by comparing 
results from two lattice volumes, 24^ x 48 and 32^ x 64, at /3 = 1.90 at one value of the sea 
quark mass {afisea = 0.0040) we notice no systematic finite volume effect on the values of 
the bag parameters. The results for Bi agree within at most one standard deviation in the 
worst case, while for the majority of the cases they are practically indistinguishable. 

3 Final results and error budget 

In this section we present the final results for the bag-parameters and we discuss the error 
budget of statistical and systematic uncertainties. 
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Figure 7: Combined chiral and continuum extrapolation for the five bag- 

parameters, Bi, renormalized in the MS scheme of [ [sT) at the scale of 3 GeV. Left and 
right panels correspond to Ml-type and M2-type four-fermion RCs, respectively, following 
the nomenclature of Ref. | [73| . In each panel open circles and stars represent the value at the 
physical point corresponding to the linear and NLO HMChPT fit, respectively. For Bi , B 2 
and B^ the polynomial (linear) and the HMChPT fits are practically indistinguishable. 


In our analysis we combine results obtained by using several possible ways to account 
for systematic effects related to the RCs determination, chiral extrapolation and discretisa¬ 
tion uncertainties. We have analysed a number of 32 estimates for Bi and 64 estimates for 
Bi with i > 1, see below for details. 

In particular, we have examined in detail the impact on the final values of the bag- 
parameters of various possible sources of systematic error related to the computation of 
the RCs. We would like to mention that a large part of the uncertainties in the RI-MOM 
calculation of the RCs affects the cutoff systematics in the error budget. 

As described in Appendix we have computed the 5x5 four-fermion RCs in the 
RI'-MOM scheme using two different methods to deal with cutoff effects, which following 
Ref. ||73| we label Ml and M2. The Ml method consists in removing 0((ap)^) effects 
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in the matrix elements used to extract the RCs by performing a fit in an appreciably large 
fixed window of the (ap)^ momentum variable. In the M2 method the RCs are determined 
as weighted averages of RCs estimators over a p^-interval (fixed in physical units) and 
common to all the gauge configuration ensembles. To control possible systematic effects 
due to the choice of the momentum interval two sets of momentum intervals have been 
compared leading to fully compatible results. 

Note that in the mixed action setup of | [5^ the off-diagonal wrong chirality mixing 
elements of the 5x5 renormalization matrix are only O(a^) cutoff effects. If the latter are 
ignored, the lattice RCs matrix shows the same mixing pattern as in the formal continuum 
theory. To check to what extent discretisation systematics can affect the final values of the 
bag-parameters, we compared the numbers obtained by simply ignoring the off-diagonal 
RC matrix elements with what one gets by including O(a^) mixing effects. 

The analysis of systematic uncertainties due to the use of polynomial and (HM)ChPT 
fit ansatze, see Eqs ( [T7| - [T9] ), is performed with reference to the so-called “golden” bag- 
combinations 


- 1 ’ 




( 20 ) 


G24 — ^^2-84, G21 — 


B2 

Bi 


Since these quantities are constructed in a way that chiral logarithmic terms cancel up to 
NLO|^ they are expected to follow an almost linear behaviour as a function of fig. Using 
the parametrization (20) we obtain additional estimates for i? 2,...,5 without having to fit 
chiral logarithmic behaviours. 

Finally, in order to estimate systematic uncertainties due to cutoff effects for i? 2,...,5 
we have also carried out the scaling analysis of quantities which are found, empirically, to 
be affected by reduced discretisation errors. Therefore, if the M 1-type RCs are employed, 
in the K case we consider 


Bi X B 2 , B 2 /B 3 , B^/Bi, Bi/B^, (21) 

while in the D case we take 

B2/B1, B2/B3, B^fB^, B^/Bg,. (22) 

If the M2-type RCs are employed, in the K case they are 

Bi X B 2 , Bi X B 3 , Bi X i?4, Bi/B^ , (23) 

while in the D case we take 

Bi X B 2 , Bi X i?3, Bi X i?4, Bi X . (24) 


Naturally in this kind of analysis all the BSM bag-parameters will eventually turn out to be 
expressed in terms of Bg which, however, among all the others is the quantity that is less 
affected by discretization effects. 

^Strictly speaking in the case of D this is not so for the combination G 24 as the NLO logarithmic terms 
do not cancel out completely. However, since the G 24 data show anyway a good linear behaviour vs. the light 
quark mass, we tried a linear fit ansatz even in this case. 
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B^/B/i, Bi X B 2 




0.00 0.01 0.02 0.03 0.04 0.00 0.01 0.02 0.03 0.04 

(li (GeV) fi£ (GeV) 

Figure 8: Combined chiral and continuum extrapolation in the case of the com¬ 

binations defined in Eq. ( |2T] ). We use Ml-type RCs in the MS scheme of | |^ at 3 GeV. For 
the combinations shown in the right panels the polynomial (linear) and the NLO ChPT fit 
ansatz coincide. 


To summarise we have carried ouj^32 kinds of analysis for Bi and 64 for the BSM 
bag-parameters, i? 2 ,..., 5 - 

Statistical errors have been evaluated using the jackknife method. We have verified 
fhaf for all fhe gauge configurafion ensembles 16 jackknife bins are enough fo have aufo- 
correlafions well under confrol. Fif cross correlafions are faken info accounf by generafing 
1000 boofslrap samples for each gauge configurafion ensemble. The RCs compufafion has 
been performed on a differenl sef of = 4 gauge configurafion ensembles. The error 
on each RC has been propagated assuming RCs fo be gaussian disfribufed wifh fhe cenfral 
values and fhe sfandard deviafions reported in Tables]^ and 

Our fofal sfafisfical uncerfainly includes fhe sfafisfical errors on fhe bare mafrix el- 
emenfs, fhe sfafisfical uncerfainly of fhe RCs and fhe propagated error coming from fhe 
combined conlinuum and chiral fif exlrapolafion. 

For each bag-parameter fhe cenfral value is determined by fhe average over fhe corre¬ 
sponding sef of resulls. Nole fhaf, since all our analyses are characterised by comparable 
fif qualify, we combine fhe resulls from differenl analyses assuming fhe same weighl for all 
of Ihem. Therefore for fhe final cenfral values as well as for fhe sfafisfical and syslemalic 

^The total number of different analyses for the SM bag-parameter Bi we have considered is given by the 
product 32 = 2x4x2x2. These numbers refer to the two fit ansatze for the chiral extrapolation, the four 
ways of combining the Ml and M2 kinds of RCs estimates needed for the renormalization of the four- and 
two-fermion operators, the two choices of the p^-interval, and finally a factor of 2 for including the off-diagonal 
scale independent O(a^) matrix elements in the construction of the renormalized operators or setting them 
equal to zero. As for the number of BSM bag-parameter analysis, owing the alternative ways of parametrizing 
chiral (Eq. |20^ and lattice spacing fit anzatze (Eqs. l|2TJ, ([2^ and Eqs. ( |22^ , \2A\), the above number must 
multiplied by 2, thus giving in total 64 kinds of analysis. 
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uncertainties we make use of the formulae (as already done in Ref. |[65|): 


X 


a 


2 


N 


N 

1 

N 




2 = 1 
N 


+ 


2=1 




(25) 

(26) 


where Xi and cjj are the central value and the variance of the i-th analysis and N is the total 
number of analyses, i.e. N = 32 for Bi and = 64 for From the first term of 

the r.h.s. of Eq. (26 1 we read off the statistical error while the second, which represents the 
spread among the results of different analyses, provides an estimate for the total systematic 
uncertainty. 

In Tables and 1^ we have collected our final resulfs for Bi {i = 1,..., 5) evaluated 
in fhe MS and RI' schemes. The final uncertainly is given by summing in quadralure Ihe 
slafislical and fhe sysfemalic errors following Eq. ( |2^ . 

In Eig. I^and Eig. we illustrate the distribution of the results for each Bi in the K 
and D case, respectively. Had we chosen the median of the results to represent the central 
value, we would have obtained numbers fully compatible (within better than one -statistical- 
standard deviation) with the results collected in Tables [T] and It has also been checked 
that the width of the interval which selects the 68% of the area around the average (or the 
median) is in all cases very close to the value provided by Eq. (261. We consider this a nice 
test of the validity and usefulness of our way of estimating the total erroij^ 

In Tables]^ and 1^ we report the detailed error budget of our determination of the K and 
D bag parameters. The numbers represent the percentage of the main sources of uncertainty 
in our calculation. The total percentage error is reported in the last rows of Tables]^ and 

Under the label “staU-fit-i-RCs” we lump together the error coming from the statistical 
uncertainties of correlators, the interpolation/extrapolation of the simulated quark masses 
to the physical values, the extrapolation to the continuum limit, as well as the statistical 
uncertainties of the RCs. 

Under the label “syst. chiral” we give our estimates of the chiral fit uncertainty. This 
has been determined using the different ways we have used to perform the chiral extrapola¬ 
tion, namely comparing results coming from the use of Eqs. ( [TT] ), ( [T^ and ( [T^ as well as 
Eqs. ( |^ . 

In the row labeled “syst. discr.” we give our final esfimafe of fhe sysfemalic uncerlain- 
lies related fo fhe choice of fhe fil ansafz of discrefisafion arlifacls. The uncertainly is laken 
as fhe spread of fhe bag parameter resulfs oblained wilh fhe use of Eqs. ( |2T] ) lo (24l and of 
fhe differenl ways (Ml or M2) of computing fhe relevanl RCs. 

Einally, in fhe Iasi row of each fable we quofe our esfimafe for fhe sysfemalic uncer- 
lainly due lo fhe perlurbafive malching of fhe RT and MS schemes. We recall lhal anoma¬ 
lous dimensions of fhe four-fermion operalors are known up lo NEO. Our associaled sys- 
lemafic error has been eslimafed by considering fhe difference belween fhe values oblained 


*The fit quality for the vast majority of the 32 analyses for Bi and the 64 analyses for Bi= 2,...,5 is good 
while only for a small number of cases - in particular, 6 (7) out of 64 analyses for B 2 and 10(11) out 64 analyses 
for B 4 for the neutral Kaon (D) case - we noticed poor fit quality. Nevertheless, we decided to attribute to the 
results of all analyses the same weight since this choice led to somewhat more conservative estimates of the 
systematic uncertainties. Infact if we had opted for a x^-weighted analysis strategy the shifts of the central 
values for all Bi would have been minimal and the finally estimated systematic errors smaller (by a few percent 
for Bi=i, 2,4 up to almost 30% for 5 ) than the ones presented in Tablesj^andj^ 
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Figure 9: Distribution of Si ,.,,,5 results for the neutral iF-mixing renormalized in the MS 


seheme of Ref. |311 at the seale of 3 GeV. The solid vertieal line marks the central value 


(average) while the gray band indicates the systematic error determined from Eq. (26 1 . 
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Figure 10: Same as in Fig.j^but for the neutral S-mixing. 


at NFO and FO at the scale of 3 GeV for each one of the bag-parameters and multiplying it 
with the value (~ 0.25) that a^®(3 GeV) takes at the same scale. 

Finite volume effects, as mentioned in the previous section, are practically negligible 
at the level of our precision. 

In Figs and we graphically show the error budget associated to the lattice com¬ 
putation i.e. without including the systematic error due to the perturbative conversion from 
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K°-K^ 

source of uncertainty (%) 


B2 

Bs 

Bi 

B5 

stat-i-fit-i-RCs 

2.5 

2.4 

2.7 

2.7 

5.4 

syst. chiral 

0.8 

1.3 

1.1 

1.8 

2.6 

syst. discr. 

2.0 

4.7 

5.8 

3.8 

4.1 

RT-MS matching 

0.5 

2.5 

1.8 

3.9 

2.3 

Total 

3.4 

6.0 

6.7 

6.3 

7.6 


Table 4: Full error budget of the estimates for the neutral iF-mixing. 


Z?'' - 

source of uncertainty (%) 

Bi 

B2 

B3 

Bi 

B5 

stat-i-fit-i-RCs 

2.9 

2.9 

4.4 

3.5 

5.1 

syst. chiral 

0.2 

0.4 

0.6 

2.3 

2.6 

syst. discr. 

2.1 

4.3 

6.7 

3.8 

3.3 

RT-MS matching 

0.5 

2.5 

1.7 

3.9 

1.1 

Total 

3.6 

5.8 

8.2 

6.9 

6.7 


Table 5: Full error budget of the estimates for the neutral D-mixing. 

RT to MS scheme, see Tables an di 



■ stat+fit+RCs 

■ syst. chiral 

■ syst. discr 


Figure 11: Graphical representation of the error budget owing to the lattice computation 
{i.e. without including the estimate for the systematic uncertainty due to the perturbative 
matching between the RT and MS schemes) for the K bag-parameters, as reported in Ta¬ 
ble Q 
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Figure 12: Same as in Fig[^for the D bag-parameters. 
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Appendix A Computational setup for the RCs 


As we use a mass-independent renormalization scheme, the calculation of the RCs of two- 
and four-fermion operators and in particular of operators with non-vanishing anomalous 
dimension must be performed in the massless quark limit. 

Eor this purpose we have produced dedicated sets of Nf = 4 Wilson twisted-mass 
degenerate dynamical quark gauge configurations with the same gluon action as the one 
used in the non-degenerate case and for a number of moderately light sea masses. Eor each 
ensemble with given sea quark mass parameters we have also computed the RCs estimators 
at several values of the valence parameters. Naturally the RCs computed with either Nf = 
2 -|- 1 -|- 1 or Nf = 4 ensembles would yield identical numbers in the chiral limit. 

The Nf = 4 ensembles are generated at values of the twist angle somewhat different 
from 7r/2 (maximal twist), i.e. at mo / rricr- The reason is that for small values of mo—mcr 
large autocorrelation times have been noticed for simulations performed at two out of the 
three values of the inverse gauge coupling (/3 = 1.90 and 1.95) we use. Although an off 


maximal-twist setup does not lead to automatic 0(a)-improvement, one can prove |54| 


that for any hadronic observable the average over results obtained at opposite values of 
the PCAC quark mass is actually 0(a)-improved. Naturally, the need of performing the 
average leads to doubling the CPU time-cost of the calculation, which however remains 
quite affordable as we are dealing with simulations at non-zero standard Wilson and twisted 
mass. 


In the Appendix A of Ref. 1651 we have presented in detail the Nf = 4 operator 
renormalization procedure for the cases of quark field and quark bilinears. Nevertheless, 
for the reader’s convenience and to fix our notations, we briefly summarise here the main 
parts of our Nf = 4 computational setup. 

We employ the Iwasaki action for the gluons while the Nf = 4 fermionic action in the 
so-called twisted basis reads 


Qsea 

‘^tm 


= a 


xj 




7 • V + lUer + {ml^j - mcr) + 


Xf 


(A.1) 
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where / = n, d, s, c, 7 • V = 7 /i(V^ + V* )/2 and Wcr = —(a/2)V* + nicr- We have 

also set 


„sea _ sea 

— ~^u ) 


^sea _ _^sea 


,fiea _ , sea _ ,,sea _ ,.sea 

Mu — l^d — ^^s — Me 


M 


(A.2) 


Note that the form of the action (A.2 1 guarantees the positivity of the fermion determinant. 
The valence fermion action takes the form 


cval _ „4 -val ^ ®V7*V7 i , .•„val,,val , 

b -a 2_^Xf 7 • V -+ mgj + zrj /ij 75 




xf- 


(A.3) 


In our notations the sea and valence sectors the various -Wilson parameters take 

values equal to ±1 and the twisted masses are non-negative quantities. 

In Table we report the simulation details and the quark mass parameters relevant for 
the Nf = 4 gauge ensembles defined above. For each value of the sea (twisted) quark mass, 
we have generated two gauge ensembles which are denoted by the letter “m” or “p” 
added to their label, and correspond to (nearly) opposite values of the PCAC quark mass, 

^-yyjSea 

ampcAc- 

Moreover, for each of the sea gauge ensembles quark propagators have been computed 
for a number of valence quark twisted masses, a/r™*. The measured value of the valence 
PCAC quark mass, for each sea ensemble of the “m” or “p” type is given in the 

last column of Tabled 

Based on Ref. the definition of the renormalized quark mass parameters in our 
partially quenched setup is 


tg 0 sea,val 


= Zp^ 

ry S0a,val 

Za ITT-pc AC 


^(Z^m™^)2 + (^-a,val)2 


sea,val sea,val ’ 

' f M 


(A.4) 


where Za is the RC of the (flavor non-singlet) axial current and rrfpc^AC denotes the PCAC 
quark mass computed from correlators in the sea and valence sector, respectively. The an¬ 
gles 6 ^^^ and are determined from the formulae cos( 0 j'^’^^™*) = 7 -®®'^Aa^^sea/vai 

and _ ZArn^pc]^’ respectively. If convenient, quark mass parame¬ 

ters in the valence sector may be chosen to be different from their sea counterparts. 

Our RC-estimators are evaluated at the values = {2n/L^) of the momenta, 
where 






([0,2], [0,2], [0,2], [0,3]) 

([2,3], [2, 3], [2, 3], [4, 7]), for/3 = 1.95, 

([0,2], [0,2], [0,2], [0,3]) 

([2,5] , [2,5] , [2,5] , [4,9]), for/3 = 1.90 and 2.10 (A.5) 


with Lfj_ the lattice size along the direction /x (with L 4 = T and Ti 2,3 = L). Quark fields 
obey anti-periodic time boundary conditions implemented by a constant shift, Ap 4 = vr/L 4 , 
in the time component of the four-momentum. Notice also that our final analysis of the RCs 
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'^^PCAC 


^sea 


cnri^^^ 

“”tcac 




13 = 1.90 (L = 24, T 

' = 48) 


A4m 

0.0080 

-0.0390(01) 

0.0285(01) 

-1.286(01) 

{0.0060, 0.0080, 0.0120, 

-0.0142(02) 

A4p 


0.0398(01) 

0.0290(01) 

-hl.291(01) 

0.0170, 0.0210 ,0.0260} 

-h0.0147(02) 

A3m 

0.0080 

-0.0358(02) 

0.0263(01) 

-1.262(02) 

{0.0060, 0.0080, 0.0120, 

-0.0152(02) 

A3p 


0.0356(02) 

0.0262(01) 

-hi.260(02) 

0.0170, 0.0210 ,0.0260} 

-hO.0147(03) 

A2m 

0.0080 

-0.0318(01) 

0.0237(01) 

-1.226(02) 

{0.0060, 0.0080, 0.0120, 

-0.0155(02) 

A2p 


-h0.0310(02) 

0.0231(01) 

-hi.218(02) 

0.0170, 0.021 0,0.0260} 

-h0.0154(02) 

Aim 

0.0080 

-0.0273(02) 

0.0207(01) 

-1.174(03) 

{0.0060, 0.0080, 0.0120, 

-0.0163(02) 

Alp 


-h0.0275(04) 

0.0209(01) 

-hi.177(05) 

0.0170, 0.021 0,0.0260} 

-h0.0159(02) 




13 = 1.95 (L = 24, T 

' = 48) 


Blm 

0.0085 

-0.0413(02) 

0.0329(01) 

-1.309(01) 

{0.0085, 0.0150, 0.0203, 

-0.0216(02) 

Blp 


-h0.0425(02) 

0.0338(01) 

-hl.317(01) 

0.0252, 0.02 98} 

-h0.0195(02) 

B7m 

0.0085 

-0.0353(01) 

0.0285(01) 

-1.268(01) 

{0.0085, 0.0150, 0.0203, 

-0.0180(02) 

B7p 


-h0.0361(01) 

0.0285(01) 

-hl.268(01) 

0.0252, 0.02 98} 

-h0.0181(01) 

B8m 

0.0020 

-0.0363(01) 

0.0280(01) 

-1.499(01) 

{0.0085, 0.0150, 0.0203, 

-0.0194(01) 

B8p 


-h0.0363(01) 

0.0274(01) 

-hl.498(01) 

0.0252, 0.02 98} 

-hO.0183(02) 

B3m 

0.0180 

-0.0160(02) 

0.0218(01) 

-0.601(06) 

{0.0060,0.0085,0.0120,0.0150, 

-0.0160(02) 

B3p 


-h0.0163(02) 

0.0219(01) 

-h0.610(06) 

0.0180,0.0203, 0.0252,0.0298} 

-h0.0162(02) 

B2m 

0.0085 

-0.0209(02) 

0.0182(01) 

-1.085(03) 

{0.0085, 0.0150, 0.0203, 

-0.0213(02) 

B2p 


-h0.0191(02) 

0.0170(02) 

-hi.046(06) 

0.0252, 0.02 98} 

-h0.0191(02) 

B4m 

0.0085 

-0.0146(02) 

0.0141(01) 

-0.923(04) 

{0.0060,0.0085,0.0120,0.0150, 

-0.0146(02) 

B4p 


-h0.0151(02) 

0.0144(01) 

-h0.940(07) 

0.0180,0.0203, 0.0252,0.0298} 

-h0.0151(02) 




l3 = 2.10 (L = 32, T 

' = 64) 


C5m 

0.0078 

-0.00821(11) 

0.0102(01) 

-0.700(07) 

{0.0048,0.0078,0.0119, 

-0.0082(01) 

C5p 


-h0.00823(08) 

0.0102(01) 

-h0.701(05) 

0.0190,0.0242 ,0.0293} 

-h0.0082(01) 

C4m 

0.0064 

-0.00682(13) 

0.0084(01) 

-0.706(09) 

{0.0039,0.0078,0.0119, 

-0.0068(01) 

C4p 


-h0.00685(12) 

0.0084(01) 

-h0.708(09) 

0.0190,0.0242 ,0.0293} 

-h0.0069(01) 

C3m 

0.0046 

-0.00585(08) 

0.0066(01) 

-0.794(07) 

{0.0025,0.0046,0.0090,0.0152, 

-0.0059(01) 

C3p 


-h0.00559(14) 

0.0064(01) 

-h0.771(13) 

0.0201,0.0249 ,0.0297} 

-h0.0056(01) 

C2m 

0.0030 

-0.00403(14) 

0.0044(01) 

-0.821(17) 

{0.0013,0.0030,0.0080,0.0143, 

-0.0040(01) 

C2p 


-h0.00421(13) 

0.0045(01) 

-h0.843(15) 

0.0195,0.0247 ,0.0298} 

-h0.0042(01) 


Table 6: Simulation details and quark mass parameters of the Nf = 4 gauge ensembles 
employed in the RCs eomputation. 


estimators has been earried out using “demoeratie” four-momentum values that satisfy the 
eondition 

~4 

A4 (p) = < 0.29, (A.6) 


with 


p^ = - sin(ap^) . 
a 


(A.7) 
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In order to come up with smoother discretisation errors we have subtracted from the Green’s 
functions entering the RI-MOM computation the perturbative cutoff effects up to order 
O(a^p^) (see Refs |74 751). 

As it has been stated above, following the general proof given in Ref. | |54| , which 
in the Appendix A.2 of Ref. 1651 has been exemplified for the case of quark bilinear RCs 
estimators, the average over RCs estimators computed on ensembles produced with opposite 
values of the sea and valence PCAC quark mass enables to remove all the odd integer power 
cutoff effects and hence the 0{a) discretisation errors ^ Based on the definition of the angle 
6 in terms of the PCAC quark mass given in Eq. (A.4 1 , we generally refer to this procedure 
as 0-average 0(o)-improvement. 

The evaluation of the RCs for the quark bilinear operators, namely Za, Zy, Zp and 
Zs as well as the RC of the quark wave function, Zq, has been done in Ref. |65|. Our 
final numbers af each value of /3 have been labeled as Ml or M2 RCs. As explained in 
defail in | |65| , fhey correspond fo differenl ways in which fhe culoff effecls are freafed. For 
convenience all fhe resulfs are again reporfed in fhe Table of Appendix of fhe presenf 
work. Nofe fhaf as for Zy, in fhe presenf analysis we have made use of fhe much more 
precise Ward Takahashi Idenfify deferminafion |731. 


Appendix B RI-MOM computation of RCs of the four-fermion 
operators 


The RT-MOM renormalizafion procedure we used for fhe four-fermion operafors has been 
explained in fhe Appendices A and B of Ref. |12|. From fhe concepfual and operafional 
poinf of view a greaf par! of fhe compufafional defails are very similar befween fhe Nf = 2 
case of Ref. 112| and fhe presenf = 4 case, excepf for fhe facf fhaf in fhe laffer one 
has fo compufe RCs estimators in fhe “m” and “p” ensembles separately before faking fheir 
0-average. In Ibis section we will fix our nofafions making extensive use of fhe descripfion 
of Ref. 112|. We will however recall some essenfial poinfs of fhe compufafion in order fo 
make easier for fhe reader fo follow fhe presenfafion of fhe analysis and our resulfs. 

In compufing RCs if is convenienf to work in a basis where fhe operafors wifh 
i = 2,..., 5 defined in Eq. ([^, are Fierz fransformed as suggesfed in Ref. m- Using 
here a generic labeling (fhaf can be obviously adapfed fo fhe case of inferesf (1, 2, 3, 4) —)■ 
{h, h', £')) fhe operator basis now read^ 


Qn±] = 2{ ([ 917 /^ 72 ][037/^74] + [ 717 ^ 7592 ][ 737 /^ 7574 ]) ± (2 o 4 )} 

Q¥[±] = 2 { ([ 917 ^ 72 ] [ 037 /^ 74 ] - [ 717 ^ 7592 ] [ 737 /^ 7574 ]) ± (2 o 4 )} 

Q ¥[±] = 2 { ([9172] [7374] - [717572] [737594]) ± (2 o 4 )} 

Q ¥[±] = 2 {([ 7 i 72 ][ 7374 ] + [717572] [737574]) ± (2 o 4 )} 

Q^it] = 2{([7ic^Mi'72][730-Mi'74]) ± (2 O 4)} (B.l) 


where color indices are mean! fo be confracfed wifhin each square parenfhesis, “MA” sfands 
for “Mixed Action” and £7^,^ = [7^,7i,]/2. The fransformafion mafrix befween fhe fwo 


^The proof for the case of the RCs of the four-fermion operators required in our mixed action setup (see 
Section!! 


10^ 


2.3 I is closely analogous to the one for the case of quark bilinear RCs. 

The quark fields qi,q 2 , q 3 and 54 are valence fields with the lattice action specified in Eq. l[^ - i.e. they are 


written in the physical basis of maximally twisted LQCD. 


22 

















operator bases (Eqs ([^ and ( |B.1| )) is given by 

/ 1 


O. 


MA _ 


*[±] 


ij ^j[±] ’ 


aM = 


0 

0 

0 


0 0 0 0 

0 0 1 0 

0 0 =f1/2 ±1/2 

0 10 0 

0 0 


(B.2) 


\ 0 ± 1/2 0 0 0 / 

To simplify the notation, in the rest of the Appendix we drop the superscript “MA” 


and subscript “±” and denote the operators (B.l I simply with the symbol Qi. Then in a 


self-evident matrix notation the renormalization pattern of the bare operators takes the 
form 

= Zq [I ± A] (B.3) 

where Zq is a scale-dependent renormalization matrix which has the same block-diagonal 
form as the formal continuum one. The wrong chirality mixings are parametrized by A 
which is a sparse off-diagonal and UV-finite matrix with the structure 


0 

Ai2 

Ai3 

Ai4 

Ai5 

A21 

0 

0 

A24 

A25 

A31 

0 

0 

CO 

<1 

A35 

A41 

A42 

A43 

0 

0 

A51 

A52 

A53 

0 

0 


(B.4) 


The renormalization pattern ( |B.3| ) can be proved (following Appendix A of Ref. 1121) for the 
renormalization of Q defined out of maximal twist, for both positive and negative mpcAC 
masses - we refer to the ensembles “p” and “m” discussed in the previous Appendix. At 
this level of course lattice artifacts are still 0(a). By the 0-averaging procedure discussed 


above, however, we obtain 0(a) improved RC-estimators for which Eq. (B.31 holds with 
only O(a^) lattice artifacts. These RC-estimators are used for renormalizing the bare matrix 
element we computed at maximal twist. 

Eor completeness we summarise the main technical points of the RC calculations. We 
start by computing in the Eandau gauge the four-point Green’s function of the operators, Qi 
between external quark states with the momenta given in Eqs (A.5 1 -( A.7| ). 

The RT-MOM renormalization condition is imposed by requiring the projected ampu¬ 
tated functions be equal to their tree-level value. This last step is conveniently and com¬ 
pactly implemented with the construction of the so-called dynamic matrix defined by fhe 
equafion D = PA, where P and A are fhe mafrices of fhe spin projectors and fhe ampufafed 


Green’s functions, respectively; see also Appendix B.l of Ref. 1121. 

The furlher steps of fhe analysis are fhe following. 

- We subfracf from fhe dynamic mafrix and fhe quark form facfor (relevanf for Zg) fhe per- 
furbafive Oja^p^^^g^) corrections computed in fhe massless fheory |74 The boosted 
coupling is defined as ^^oost ~ 6/(/3(P)) where {P) is fhe average plaqueffe value. 

- We subfracf from fhe projecfed ampufafed four-fermion correlafors fhe confribufion of 


fhe Goldsfone boson pole (GB-pole) af each value of fhe momenfa defined in Eqs (A.51 


(A.7 1 . This sfep is performed by carrying ouf fhe chiral limif exfrapolafion in fhe va- 
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lence sector. The GB-pole contribution is realised by the presence of terms that go likep^ 
(l/TTip^)", being each one of them suppressed by corresponding powers of the inverse 
square momentum, (1/p^)”. For each sea quark ensemble and for each value of the mo¬ 
mentum we fit each one of the dynamic matrix elements in terms of the pseudoscalar 
mass in the valence sector, using the ansatz 

Ai(aV;"ipJ = + ^(aV)"iL +C(aV)/"rps + (B.5) 


In Fig. 13 we show the GB-pole subtraction of the block-diagonal elements of the dy¬ 
namic matrix corresponding to the continuum-like matrix elements in the case of the 
coarsest lattice (/3 = 1.90) at relatively small momenta where systematic effects are ex¬ 
pected to be larger. We observe, as it is also expected, that the GB-pole contribution is 
getting suppressed as the value of the momentum increases. This feature can be noticed 
comparing Figs [T^ and both referring to /3 = 1.90 at two different values of the mo¬ 
mentum. Moreover typical GB-pole subtraction fits in the finest lattice are given in the 


panels of Fig. 15 


Important information concerning the double GB-pole subtraction can be revealed from 
plots as for example those presented in the panels of Fig. 16 which refer to the ensem¬ 
bles Alp and Aim of the coarsest lattice spacing (/3 = 1.90). We form the product of 
with the fit parameter, 'D{a?p^), of the double GB-pole term {cf. Eq. (B.5 l), then 
plotting it against Two observations are in order. First, we notice that the product 

V{a^p^) X takes almost constant value for large enough values of momentum 

that lie in the momentum intervals we have used in order to extract our RCs estimates. 
This finding also serves as a confirmation of the good fit quality in analysing the double 
GB-pole term. Second, we find that the double GB-pole contribution is negligible for the 
case of Dll, in nice agreement with theoretical expectations, see Ref. | [6^ , whereas it is 
different from zero for several of the Dij with i,j > 1. So in our final analysis we have 
adopted single GB-pole subtraction (i.e. set P = 0 in the fit ansatz of Eq. ( |B.5| )) for the 
former and single and double GB-pole subtraction (cf. Eq. (B.51) for the latter cases. 

After applying the 0-average of the RC-estimators at each /3 and at each value of mo¬ 
mentum defined by Eqs (A.6l and (A.7i, which is required in order to achieve 0(a)- 
improvement, we carry out the sea chiral limit for each element of the 5x5 renormaliza¬ 
tion matrix. A simple polynomial (linear) fit ansatz fits the data smoothly. In Eig.fTVjwe 


illustrate an example of the chiral extrapolation in the sea using almost the same value of 
momentum at the three /?’s. 

Having performed the chiral limit extrapolations our estimates for the RCs take the form 
(P^'i {P^'i = 0), where the first argument refers to 

the scale. Moreover we obtain the scale independent off-diagonal elements, namely 
{{ap)^) = A^^'((ap)^, = 0). The estimators (ap)^) can be 


evolved to a common scale po using the running formula for the operators Qi known 
up to NEO 131 77 1 obtaining thus estimates of the form (pq; {ap)‘^). To be able to 
carry out a controlled study of the systematic discretisation errors on the renormalized 
bag-parameters, we apply the two methods proposed in Ref. | 73 |. We recall that each 
of these two methods, called Ml and M2, prescribes a different treatment of the cutoff 
effects. Method Ml consists in fitting (pg; {ap)'^) to the linear ansatz 


io,pY) = + Xij X (ap) 


(B.6) 


"For each combination of two valence quark masses in each gauge ensemble of Tablewe have computed 
and used the corresponding value of the pseudoscalar mass, nips- 
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Figure 13: Fitting procedure for the GB-pole subtraction on the block diagonal elements 
of the dynamic matrix. We show two examples from the coarsest lattice, in particular for 
the ensemble A2m (left) and A2p (right), at a relatively small value of momentum, namely 
(ai 


2 ~ 1.57. 
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Figure 14: Same as in Fig. 
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13 


but at a larger value of momentum, namely (ap)^ ~ 2.19. 
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Figure 15: Same as in Figj^at /3 = 2.10, for the ensemble C2m (left) and C2p (right) and 
(ap)^ ~ 1.57. 
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V{a^p^) X {a^p^)^ P(a^p^) x {a^p^)^ P(a^p^) x {a^p^) 



Figure 16: For each of block-diagonal matrix element we form the product of {a^p^)'^ with 
the fit parameter 'D{a^p^) that is associated to the term of the double GB-pole {cf. Eq. (B.51) 
and we plot it against (a^p^) for the ensembles Alp and Aim of /3 = 1.90. 


/3 = 1.90 


P = 1.95 


P = 2.10 



x 10“ 


(aA/sea)2 


^ 

x 10“ 


Figure 17: Example of the sea chiral extrapolation of the RC estimators at three values of 
/3 computed at (ap)^ ~ 1.57. 


in some large momentum region. Eor better controlling the systematics we have made 


two choices, namely {app G [1.5, 2.2] and {app G 


subtraction of the perturbative 0{a‘^g‘^ '' 


[1.8, 2.2]. As expected, thanks to 
show very smooth depen- 


effects, the slopes Ajj 

dence on /3. In fact we parametrize the slopes in Eq. Jb.6i as Ajj = a|°^ + A-j^p^ and 
we perform a simultaneous linear extrapolation to {app = 0 at all values of j3. Then we 
may convert the extrapolated results (po) to any scale and scheme as for example 


the MS scheme using NEO running. In Eig. 18 we show the best linear fits of Zij (in 


the MS scheme at the scale of 3 GeV) using Eq. (B.6l. M2 method works in a quite 
different way from M1. It consists in getting the average RC value in a narrow window 
of momenta which is fixed in physical unifs for all values of /3. We have carried ouf fhe 
M2-lype analysis for fwo choices of momenfum interval, namely p^ G [10 : 13] GeV^ 
and p^ G [11 : 14] GeV^. These rafher high values of momenfum offer fhe possibilify fo 
gain more confidence in fhe absence of hadronic sfafe confaminafions and in fhe validify 
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of the RC-evolution using the NLO anomalous dimension at the price, however, of taking 
no special care to end up with reduced O(a^) cutoff effects. 



Figure 18: Best linear fits {cf. Eq. B .61 in the momentum interval (ap)^ G [1.8 : 2.2] for the 
block-diagonal RCs, Zij, at each value of /?, renormalized in the MS scheme at the scale of 
3GeV. 


In order to get an immediate view of the impact of the RCs cutoff effects on the contin¬ 
uum limit estimates for the bag-parameters, we have performed a scaling test for all Bi 
(f = 1,..., 5) which, for the case of the neutral iF-mixing, is illustrated in Figs[T^ By 
using Ml or M2-type RCs we get the bag-parameter estimates at some fixed reference 
value of fhe lighf quark mass and fhen we exfrapolafe fhem linearly in fo fhe confinuum 
limif. From fhe relevanf panels of Fig. 19 if can be noficed fhaf for bofh RCs-fypesp^fhe 


-scaling behaviour is good and fhe exfrapolafed continuum limif resulfs are compafible 
wifhin 1 or 2 sfandard deviations, depending on fhe case. 

The sifuafion presenfed in Figs[^has fo be considered as purely indicafive buf represen- 
fafive of fhe facf fhaf af some arbrifrary -reference- value of fhe lighf quark mass, using all 
four possible combinations of fhe RCs, fhe resulfs for converge fo confinuum 

limif values fhaf are compafible wifhin each ofher. We find fhis resulf very reassuring 
since Ml- and M2-lype RCs are computed in such a way fhaf fhe corresponding cufoff 
effecfs are much differenf, fhough of O(a^). In fhis sense, by using fwo fypes of RCs we 
gain confidence fhaf sysfemafic effecfs due fo fhe RCs RI-MOM compufafion and discer- 
fisafion effecfs are under confrol. We also recall fhaf for our final resulfs we do nof rely on 


'^Notice that we have considered four cases corresponding to all possible combinations of taking Ml- and 
M2-type RCs for the four- and two-fermion operators. 
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plots like the ones presented in Fig. 19, but we perform combined chiral and continuum 
limit fits as those descibed by the fit ansatze of Eqs (17)-(19) of Section 2. 

Finally, in Fig. we depict the behaviour of the scale independent off-diagonal matrix 
elements Ajj for /3 = 2.10, while in Fig. 21 we plot the final estimates for Ajj against 
a^. We note that in all cases A^ get small values around zero. For reader’s convenience 
we collect in the last row of Fig. 1^ the Aij for which we observe relatively larger values. 
As we have anticipated in Section ^ in order to take into account possible residual cut¬ 
off effects, in our final set of analyses we have included results for the bag-parameters 
computed both with = 0 and Aij ^ 0. 


Bi B2 B3 B4 B5 




Figure 19: estimates computed at the physical strange quark value and at a fixed 

reference quark mass = 12.0 MeV plotted against for the three values of the lattice 
spacing. We compare the scaling behaviour of the bag-parameter estimates computed with 
the four possible combinations of Ml- and M2-type for the four- and two-fermion RCs. 
The best linear fit in and the corresponding CL value for each RCs combination is also 
shown. 


We collect for convenience the RCs for the bilinear quark operators, calculated in 
Ref. | [M1 , in Table while in Tables and we summarise the RCs values for the four- 
fermion operators. 
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Figure 20: The behaviour of the off-diagonal mixing eoeffieients as a function of (ap)^ 
for ^ = 2.10. 
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Figure 21: computed with the M2-method against a^. 
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1.4 1.6 1.8 2.0 2.2 2.4 


{ap? (ap)^ (ap)^ (ap)^ 
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/3 

Method 

Zy 

Za 

Zp 




Ml 

0.587(04) 

0.731(08) 

0.587(08) 

0.830(14) 

0.705(05) 

1.90 

M2 

0.608(03) 

0.703(02) 

0.637(06) 

0.974(04) 

0.720(02) 


WTI 

0.5920(04) 

- 

- 

- 

- 


Ml 

0.603(03) 

0.737(05) 

0.566(05) 

0.812(09) 

0.719(04) 

1.95 

M2 

0.614(02) 

0.714(02) 

0.606(03) 

0.913(03) 

0.727(01) 


WTI 

0.6095(03) 

- 

- 

- 

- 


Ml 

0.655(03) 

0.762(04) 

0.572(02) 

0.777(06) 

0.759(04) 

2.10 

M2 

0.657(02) 

0.752(02) 

0.605(02) 

0.832(04) 

0.760(02) 


WTI 

0.6531(02) 

- 

- 

- 

- 


Table 7: Bilinear RCs published in Ref. | [^ . The scale independent Zy, Za and the scale 
dependent Zp, Zg and Zq are obtained with the methods Ml and M2. The scale dependent 
RCs are expressed in the MS scheme at the scale of 3 GeV. Zy is also obtained performing 
a very accurate computation employing the Ward-Takahashi identity (WTI), for details see 
Section 2.3 of Ref. |[7^). 


RI' (3 GeV) 

/3 = 

1.90 

/3 = 

1.95 

/3 = 

2.10 

Zij 

Ml 

M2 

Ml 

M2 

Ml 

M2 

Zix 

0.373(07) 

0.383(03) 

0.395(05) 

0.398(05) 

0.454(04) 

0.455(03) 

Z 22 

0.450(07) 

0.450(04) 

0.474(06) 

0.469(05) 

0.536(07) 

0.528(04) 

Z‘2.1, 

0.200(07) 

0.073(03) 

0.222(05) 

0.122(05) 

0.236(03) 

0.175(02) 

Zz2 

0.015(02) 

0.009(01) 

0.015(01) 

0.010(01) 

0.015(00) 

0.013(00) 

Z-iJ. 

0.247(11) 

0.382(04) 

0.237(07) 

0.337(05) 

0.285(06) 

0.335(02) 

Z44 

0.277(08) 

0.368(03) 

0.277(06) 

0.344(04) 

0.333(05) 

0.362(02) 

■^45 

-0.012(01) 

-0.010(00) 

-0.012(01) 

-0.010(01) 

-0.012(01) 

-0.011(00) 

■^54 

-0.146(05) 

-0.054(02) 

-0.166(04) 

-0.090(03) 

-0.187(03) 

-0.141(02) 

■^55 

0.435(07) 

0.403(03) 

0.471(05) 

0.439(05) 

0.552(07) 

0.533(04) 


Table 8: Typical values for the four-fermion operator RCs at three values of /?. For Ml 
method linear extrapolation to (ap)^ has been performed using data in the interval (ap)^ G 
[1.8, 2.2], while for method M2 we have used data from the narrow momentum window 
determined by p^ G [11 : 14] GeV^. RCs are expressed in the in the RT scheme at the scale 
of 3 GeV. 
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MS (3 GeV) 

13 = 

1.90 

[3 = 

1.95 

13 = 

2.10 

Zij 

Ml 

M2 

Ml 

M2 

Ml 

M2 

Zxx 

0.379(07) 

0.389(03) 

0.402(05) 

0.404(06) 

0.462(04) 

0.462(03) 

Z 22 

0.440(07) 

0.440(03) 

0.463(06) 

0.458(05) 

0.524(07) 

0.516(04) 

■^23 

0.182(08) 

0.050(03) 

0.204(05) 

0.101(05) 

0.216(03) 

0.153(02) 

■^32 

0.020(02) 

0.013(01) 

0.020(02) 

0.014(01) 

0.021(01) 

0.017(00) 

■^33 

0.293(13) 

0.453(04) 

0.281(08) 

0.399(06) 

0.339(07) 

0.398(03) 

Z44 

0.304(09) 

0.405(03) 

0.303(07) 

0.378(05) 

0.364(06) 

0.397(03) 

■^45 

-0.006(01) 

-0.004(00) 

-0.006(01) 

-0.004(01) 

-0.005(01) 

-0.003(00) 

■^54 

-0.143(06) 

-0.042(02) 

-0.163(04) 

-0.081(03) 

-0.183(03) 

-0.134(02) 

■^55 

0.460(08) 

0.426(04) 

0.497(06) 

0.464(06) 

0.584(08) 

0.564(05) 


Table 9: Same as in Table 8 but in the MS scheme of 1311 at the scale of 3 GeV. 
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